Analyze transcript data of TIER-Seq data for Synechocystis RNase E
Ute Hoffmann
AG Annegret Wilde, Institute of Biology IIIJuli 09, 2021
- 1 Aim of the analysis
- 2 DESeq2 analysis
- 3 Exploratory Data Analysis
- 3.1 Create table how many RNA features of a certain type are affected by differential expression
- 3.2 Functional Enrichment
- 3.3 Gene Set Enrichment Analyses
- 3.4 Plots of GSEA
- 3.5 GSEA for RNA features
- 3.6 GSEA for base means
- 3.7 Without taking width into account
- 3.8 When taking width of features into account
- Session info
1 Aim of the analysis
Analyze transcript level data from TIER-Seq comparison and do some exploratory data analysis.
2 DESeq2 analysis
Read in data: htseq-count files created with usegalaxy.eu
htseq_CDS <- read.delim("input/Galaxy129-[Column_Join_on_data_115,_data_123,_and_others].tabular", header=TRUE, row.names=1)
htseq_TUs <- read.delim("input/Galaxy158-[Column_Join_on_data_146,_data_147,_and_others].tabular", header=TRUE, row.names=1)
zuordnung <- read.delim("input/20210114_zuordnungIn_20210114_transcript-history.csv", header=TRUE)
row.names(zuordnung) <- names(htseq_CDS)[c(5,1,6,2,4,3,10,7,11,8,12,9)]
names(htseq_CDS) <- zuordnung[names(htseq_CDS),]$name
names(htseq_TUs) <- zuordnung[names(htseq_TUs),]$name
coldata <- read.csv("input/20210125_colData.csv", row.names=1)Do actual analysis
# create DESeq2 data object
ddsMat_CDS <- DESeqDataSetFromMatrix(countData = htseq_CDS,
colData = coldata[names(htseq_CDS),],
design = ~ condition)## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factorsddsMat_TUs <- DESeqDataSetFromMatrix(countData = htseq_TUs,
colData = coldata[names(htseq_TUs),],
design = ~ condition)## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors# run DESeq
ddsMat_CDS <- DESeq(ddsMat_CDS) ## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingddsMat_TUs <- DESeq(ddsMat_TUs) ## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingwrite.csv(data.frame(factor=ddsMat_CDS$sizeFactor), file="output/transcript_sizeFactors.csv")
write.csv(data.frame(counts(ddsMat_CDS, normalized=TRUE)), file="output/Transcript_CDS_normalizedCounts.csv")2.1 Diagnostic Plots
plotDispEsts(ddsMat_CDS, main="CDS comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")
plotDispEsts(ddsMat_TUs, main="TU comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")
pdf(file="output/DESeq2_Plots/CDS/ddsMat_CDS_DispEsts.pdf", width=4.5, height=4.5)
plotDispEsts(ddsMat_CDS, main="CDS comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")
dev.off()## png
## 2pdf(file="output/DESeq2_Plots/TU/ddsMat_TU_DispEsts.pdf", width=4.5, height=4.5)
plotDispEsts(ddsMat_TUs, main="TU comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")
dev.off()## png
## 2p <- PCA_plot(ddsMat_CDS, "RNA Features")
p 
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS_PCA.pdf", plot=p, width=9, height=9, units="cm")
p <- PCA_plot(ddsMat_TUs, "TU")
p
ggsave("output/DESeq2_Plots/TU/ddsMat_TU_PCA.pdf", plot=p, width=9, height=9, units="cm")
p <- heatmap_plot(ddsMat_CDS, "RNA Features")
p
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS_heatMap.pdf", plot=p, width=15, height=12, units="cm")
p <- heatmap_plot(ddsMat_TUs, "TU")
p
ggsave("output/DESeq2_Plots/TU/ddsMat_TU_heatMap.pdf", plot=p, width=15, height=12, units="cm")2.2 Extract Results
Extract results and use changeAnnotation_DESeq2.py and annotation_locusTags_stand13012021.csv to add annotation to tables.
# extract results
CDS_result_dWT_dIF_1h <- results(ddsMat_CDS, contrast=c("condition", "dWT_1h", "dIF_1h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dWT_dIF_1h[order(CDS_result_dWT_dIF_1h$padj),])), file="output/DESeq2_resultsTables/results_CDS-1h.tsv")
CDS_result_dWT_dIF_0h <- results(ddsMat_CDS, contrast=c("condition", "dWT_0h", "dIF_0h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dWT_dIF_0h[order(CDS_result_dWT_dIF_0h$padj),])), file="output/DESeq2_resultsTables/results_CDS-0h.tsv")
CDS_result_dWT <- results(ddsMat_CDS, contrast=c("condition", "dWT_1h", "dWT_0h")) # dWT 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dWT[order(CDS_result_dWT$padj),])), file="output/DESeq2_resultsTables/results_CDS_dWT.tsv")
CDS_result_dIF <- results(ddsMat_CDS, contrast=c("condition", "dIF_1h", "dIF_0h")) # dIF 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(CDS_result_dIF[order(CDS_result_dIF$padj),])), file="output/DESeq2_resultsTables/results_CDS_dIF.tsv")
TUs_result_dWT_dIF_1h <- results(ddsMat_TUs, contrast=c("condition", "dWT_1h", "dIF_1h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dWT_dIF_1h[order(TUs_result_dWT_dIF_1h$padj),])), file="output/DESeq2_resultsTables/results_TUs-1h.tsv")
TUs_result_dWT_dIF_0h <- results(ddsMat_TUs, contrast=c("condition", "dWT_0h", "dIF_0h")) # dWT/dIF -> higher in dWT: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dWT_dIF_0h[order(TUs_result_dWT_dIF_0h$padj),])), file="output/DESeq2_resultsTables/results_TUs-0h.tsv")
TUs_result_dWT <- results(ddsMat_TUs, contrast=c("condition", "dWT_1h", "dWT_0h")) # dWT 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dWT[order(TUs_result_dWT$padj),])), file="output/DESeq2_resultsTables/results_TUs_dWT.tsv")
TUs_result_dIF <- results(ddsMat_TUs, contrast=c("condition", "dIF_1h", "dIF_0h")) # dIF 1h/0h -> higher in 1h: higher log2FC
write_tsv(rownames_to_column(as.data.frame(TUs_result_dIF[order(TUs_result_dIF$padj),])), file="output/DESeq2_resultsTables/results_TUs_dIF.tsv")Execute in directory “DESeq2_resultsTables”: python changeAnnotation_DESeq2.py results_CDS-1h.tsv results_CDS-1h_annotated.tsv python changeAnnotation_DESeq2.py results_CDS-0h.tsv results_CDS-0h_annotated.tsv python changeAnnotation_DESeq2.py results_CDS_dWT.tsv results_CDS-dWT_annotated.tsv python changeAnnotation_DESeq2.py results_CDS_dIF.tsv results_CDS-dIF_annotated.tsv
python TUs_add_info.py results_TUs-1h.tsv results_TUs-1h_annotated.tsv python TUs_add_info.py results_TUs-0h.tsv results_TUs-0h_annotated.tsv python TUs_add_info.py results_TUs_dWT.tsv results_TUs-dWT_annotated.tsv python TUs_add_info.py results_TUs_dIF.tsv results_TUs-dIF_annotated.tsv
2.3 p-Value plots for different comparisons
pvaluePlot(CDS_result_dWT_dIF_0h, "CDS 0h")
pvaluePlot(CDS_result_dWT_dIF_1h, "CDS 1h")
pvaluePlot(CDS_result_dWT, "CDS dWT")
pvaluePlot(CDS_result_dIF, "CDS dIF")
pvaluePlot(TUs_result_dWT_dIF_0h, "TUs 0h")
pvaluePlot(TUs_result_dWT_dIF_1h, "TUs 1h")
pvaluePlot(TUs_result_dWT, "TUs dWT")
pvaluePlot(TUs_result_dIF, "TUs dIF")
## png
## 2## png
## 2## png
## 2## png
## 2## png
## 2## png
## 2## png
## 2## png
## 22.4 Filter for non-NA values
CDS_result_dWT_dIF_0h <- subset(CDS_result_dWT_dIF_0h, !is.na(CDS_result_dWT_dIF_0h$padj))
nrow(CDS_result_dWT_dIF_0h)## [1] 4318CDS_result_dWT_dIF_1h <- subset(CDS_result_dWT_dIF_1h, !is.na(CDS_result_dWT_dIF_1h$padj))
nrow(CDS_result_dWT_dIF_1h)## [1] 5293CDS_result_dWT <- subset(CDS_result_dWT, !is.na(CDS_result_dWT$padj))
nrow(CDS_result_dWT)## [1] 4968CDS_result_dIF <- subset(CDS_result_dIF, !is.na(CDS_result_dIF$padj))
nrow(CDS_result_dIF)## [1] 5510TUs_result_dWT_dIF_0h <- subset(TUs_result_dWT_dIF_0h, !is.na(TUs_result_dWT_dIF_0h$padj))
nrow(TUs_result_dWT_dIF_0h)## [1] 3620TUs_result_dWT_dIF_1h <- subset(TUs_result_dWT_dIF_1h, !is.na(TUs_result_dWT_dIF_1h$padj))
nrow(TUs_result_dWT_dIF_1h)## [1] 4009TUs_result_dWT <- subset(TUs_result_dWT, !is.na(TUs_result_dWT$padj))
nrow(TUs_result_dWT)## [1] 3853TUs_result_dIF <- subset(TUs_result_dIF, !is.na(TUs_result_dIF$padj))
nrow(TUs_result_dIF)## [1] 40092.5 Create MA and Volcano plots for different comparisons, count features up-, downregulated
2.5.1 CDS
2.5.1.1 CDS 0h
count_up_down(CDS_result_dWT_dIF_0h, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 20"
## [1] "number of features up: 94"p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dWT_dIF_0h), foldchange=0.8, padjusted=0.05, text=TRUE)
p## Warning: Removed 4204 rows containing missing values (geom_text_repel).## Warning: ggrepel: 50 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-0h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")## Warning: Removed 4204 rows containing missing values (geom_text_repel).## Warning: ggrepel: 90 unlabeled data points (too many overlaps). Consider
## increasing max.overlapsp <- MAplot_ggplot(CDS_result_dWT_dIF_0h, foldchange=0.8, y_axis_label = "Log2 fold-change(rne(WT)/rne(Ts))") + ylim(-5.5,+5.5)
p
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-0h_MAplot.pdf",plot=p, width=15, height=12, units="cm")2.5.1.2 CDS 1h
count_up_down(CDS_result_dWT_dIF_1h, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 433"
## [1] "number of features up: 443"p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dWT_dIF_1h), foldchange=0.8, padjusted=0.05, text=TRUE)
p## Warning: Removed 4417 rows containing missing values (geom_text_repel).## Warning: ggrepel: 817 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-1h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")## Warning: Removed 4417 rows containing missing values (geom_text_repel).## Warning: ggrepel: 863 unlabeled data points (too many overlaps). Consider
## increasing max.overlapsp <- MAplot_ggplot(CDS_result_dWT_dIF_1h, foldchange=0.8, y_axis_label = "Log2 fold-change(rne(WT)/rne(Ts))") + ylim(-5.5,+5.5)
p
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-1h_MAplot.pdf",plot=p, width=15, height=12, units="cm")2.5.1.3 CDS dWT
count_up_down(CDS_result_dWT, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 360"
## [1] "number of features up: 302"p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dWT), foldchange=0.8, padjusted=0.05, text=TRUE)
p## Warning: Removed 4306 rows containing missing values (geom_text_repel).## Warning: ggrepel: 618 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dWT_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")## Warning: Removed 4306 rows containing missing values (geom_text_repel).## Warning: ggrepel: 645 unlabeled data points (too many overlaps). Consider
## increasing max.overlapsp <- MAplot_ggplot(CDS_result_dWT, foldchange=0.8, y_axis_label = "Log2 fold-change(rne(WT) 1h / 0h)") + ylim(-5.5,+5.5) + scale_colour_manual(values=c("DOWN"="black", "UP"="black", "NO"="#d3d3d3ff"))## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.p## Warning: Removed 1 rows containing missing values (geom_point).
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dWT_MAplot.pdf",plot=p, width=15, height=12, units="cm")## Warning: Removed 1 rows containing missing values (geom_point).2.5.1.4 CDS dIF
count_up_down(CDS_result_dIF, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 636"
## [1] "number of features up: 686"p <- volcanoPlot_ggplot(as.data.frame(CDS_result_dIF), foldchange=0.8, padjusted=0.05, text=TRUE)
p## Warning: Removed 4188 rows containing missing values (geom_text_repel).## Warning: ggrepel: 1277 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dIF_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")## Warning: Removed 4188 rows containing missing values (geom_text_repel).## Warning: ggrepel: 1311 unlabeled data points (too many overlaps). Consider
## increasing max.overlapsp <- MAplot_ggplot(CDS_result_dIF, foldchange=0.8, y_axis_label = "Log2 fold-change(Rne-Ts 1h / 0h)") #+ ylim(-5,+5)
p
ggsave("output/DESeq2_Plots/CDS/ddsMat_CDS-dIF_MAplot.pdf",plot=p, width=15, height=12, units="cm")2.5.2 TUs
2.5.2.1 TUs 0h
count_up_down(TUs_result_dWT_dIF_0h, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 8"
## [1] "number of features up: 78"p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dWT_dIF_0h), foldchange=0.8, padjusted=0.05)
p
ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-0h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
DESeq2::plotMA(TUs_result_dWT_dIF_0h)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-0h_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dWT_dIF_0h, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()## png
## 22.5.2.2 TUs 1h
count_up_down(TUs_result_dWT_dIF_1h, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 327"
## [1] "number of features up: 309"p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dWT_dIF_1h), foldchange=0.8, padjusted=0.05)
p
ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-1h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
DESeq2::plotMA(TUs_result_dWT_dIF_1h)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-1h_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dWT_dIF_1h, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()## png
## 22.5.2.3 TUs dWT
count_up_down(TUs_result_dWT, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 268"
## [1] "number of features up: 227"p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dWT), foldchange=0.8, padjusted=0.05)
p
ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-dWT_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
DESeq2::plotMA(TUs_result_dWT)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-dWT_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dWT, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()## png
## 22.5.2.4 TUs dIF
count_up_down(TUs_result_dIF, foldchange=0.8, padjusted=0.05)## [1] "number of features down: 492"
## [1] "number of features up: 524"p <- volcanoPlot_ggplot(as.data.frame(TUs_result_dIF), foldchange=0.8, padjusted=0.05)
p
ggsave("output/DESeq2_Plots/TU/ddsMat_TUs-dIF_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
DESeq2::plotMA(TUs_result_dIF)
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
pdf(file="output/DESeq2_Plots/TU/ddsMat_TUs-dIF_MAplot.pdf", width=4.5, height=4.5)
DESeq2::plotMA(TUs_result_dIF, xlab="Mean of Normalized Counts", ylab="Log2FC(rne-OE/rne-ts)")
abline(h=c(-0.8,+0.8), col="dodgerblue", lwd=2)
dev.off()## png
## 23 Exploratory Data Analysis
3.1 Create table how many RNA features of a certain type are affected by differential expression
features <- rtracklayer::import("input/20210217_syne_onlyUnique_withFeat.gff3")
CDS_features <- subset(features, features$type=="CDS")
TUs <- rtracklayer::import("input/Kopf_4091_TUs_combined.gff3")
types <- rep(c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc"),2)
updown <- c(rep("up",9), rep("down",9))
df_features <- data.frame(cbind(type=types, updown=updown))
types <- c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc")
for(t in types){
t_feat <- subset(features, features$type==t)
subset_CDS <- subset(CDS_result_dWT, row.names(CDS_result_dWT) %in% t_feat$locus_tag)
index_up <- which(df_features$type==t & df_features$updown=="up")
index_down <- which(df_features$type==t & df_features$updown=="down")
df_features[index_up,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange>0.8)) # count number of features affected
df_features[index_down,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange<(-0.8))) # count number of features affected
df_features[df_features$type==t,"number_total"] <- length(t_feat)
}
df_features## type updown number_feat_overlap number_total
## 1 CDS up 193 3675
## 2 5UTR up 25 979
## 3 3UTR up 3 29
## 4 tRNA up 0 42
## 5 rRNA up 1 6
## 6 ncRNA up 28 318
## 7 asRNA up 52 1071
## 8 CRISPR up 0 3
## 9 misc up 0 11
## 10 CDS down 292 3675
## 11 5UTR down 26 979
## 12 3UTR down 2 29
## 13 tRNA down 5 42
## 14 rRNA down 0 6
## 15 ncRNA down 10 318
## 16 asRNA down 24 1071
## 17 CRISPR down 0 3
## 18 misc down 1 11write.csv(df_features, file="output/RNAfeatures_upDown_dWT_0h1hheat.csv")# prepare data.frame for barplot
types <- rep(c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc"),2)
updown <- c(rep("up",9), rep("down",9))
df_features <- data.frame(cbind(type=types, updown=updown))
types <- c("CDS","5UTR", "3UTR", "tRNA", "rRNA", "ncRNA", "asRNA", "CRISPR", "misc")
for(t in types){
t_feat <- subset(features, features$type==t)
subset_CDS <- subset(CDS_result_dWT_dIF_1h, row.names(CDS_result_dWT_dIF_1h) %in% t_feat$locus_tag)
index_up <- which(df_features$type==t & df_features$updown=="up")
index_down <- which(df_features$type==t & df_features$updown=="down")
df_features[index_up,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange>0.8)) # count number of features affected
df_features[index_down,"number_feat_overlap"] <- nrow(subset(subset_CDS, subset_CDS$padj<0.05 & subset_CDS$log2FoldChange<(-0.8))) # count number of features affected
df_features[df_features$type==t,"number_total"] <- length(t_feat)
}
df_features## type updown number_feat_overlap number_total
## 1 CDS up 333 3675
## 2 5UTR up 29 979
## 3 3UTR up 2 29
## 4 tRNA up 17 42
## 5 rRNA up 0 6
## 6 ncRNA up 28 318
## 7 asRNA up 31 1071
## 8 CRISPR up 3 3
## 9 misc up 0 11
## 10 CDS down 226 3675
## 11 5UTR down 41 979
## 12 3UTR down 1 29
## 13 tRNA down 0 42
## 14 rRNA down 0 6
## 15 ncRNA down 24 318
## 16 asRNA down 140 1071
## 17 CRISPR down 0 3
## 18 misc down 1 11write.csv(df_features, file="output/RNAfeatures_upDown_dWT_dIF_1hheat.csv")3.2 Functional Enrichment
go_functional_enrichment(CDS_result_dWT, write=TRUE, path_up="output/enrichment/go_enrichment/dWT_0h_1h_up.csv", path_down="output/enrichment/go_enrichment/dWT_0h_1h_down.csv")## ID Description GeneRatio BgRatio pvalue
## GO:0006310 GO:0006310 DNA recombination 7/130 27/2459 0.0003587321
## p.adjust qvalue
## GO:0006310 0.03694941 0.03694941
## geneID Count
## GO:0006310 sll0832/sll6059/sll6109/slr0181/slr5010/slr7005/slr8029 7
## ID Description GeneRatio
## GO:0008137 GO:0008137 NADH dehydrogenase (ubiquinone) activity 13/240
## GO:0048038 GO:0048038 quinone binding 13/240
## GO:0042773 GO:0042773 ATP synthesis coupled electron transport 9/240
## GO:0031470 GO:0031470 <NA> 8/240
## GO:0015977 GO:0015977 carbon fixation 8/240
## GO:0015986 GO:0015986 ATP synthesis coupled proton transport 6/240
## BgRatio pvalue p.adjust qvalue
## GO:0008137 19/2459 8.545747e-10 1.076764e-07 9.805120e-08
## GO:0048038 21/2459 5.353053e-09 3.372423e-07 3.070962e-07
## GO:0042773 14/2459 9.008581e-07 3.783604e-05 3.445387e-05
## GO:0031470 12/2459 2.579196e-06 8.124469e-05 7.398221e-05
## GO:0015977 13/2459 6.143525e-06 1.548168e-04 1.409777e-04
## GO:0015986 10/2459 1.220117e-04 2.562246e-03 2.333207e-03
## geneID
## GO:0008137 sll0223/sll0520/sll0521/sll0522/sll1732/sll1733/slr0261/slr1279/slr1280/slr1281/slr1291/slr2007/slr2009
## GO:0048038 sll0223/sll0519/sll0520/sll0521/sll0522/sll1732/slr0261/slr1279/slr1280/slr1281/slr1623/slr2007/slr2009
## GO:0042773 sll0223/sll0522/sll1732/sll1733/slr1291/slr2006/slr2007/slr2008/slr2009
## GO:0031470 sll1028/sll1029/sll1030/sll1031/sll1032/slr0009/slr0012/slr0169
## GO:0015977 sll1028/sll1029/sll1030/sll1031/sll1032/slr0009/slr0012/slr0169
## GO:0015986 sll1324/sll1325/sll1326/sll1327/slr1330/ssl2615
## Count
## GO:0008137 13
## GO:0048038 13
## GO:0042773 9
## GO:0031470 8
## GO:0015977 8
## GO:0015986 6go_functional_enrichment(CDS_result_dIF, write=TRUE, path_up="output/enrichment/go_enrichment/dIF_0h_1h_up.csv", path_down="output/enrichment/go_enrichment/dIF_0h_1h_down.csv")## ID Description GeneRatio BgRatio
## GO:0005887 GO:0005887 integral component of plasma membrane 20/280 75/2476
## GO:0004673 GO:0004673 protein histidine kinase activity 13/280 43/2476
## GO:0000155 GO:0000155 phosphorelay sensor kinase activity 15/280 54/2476
## pvalue p.adjust qvalue
## GO:0005887 0.0001447987 0.02157500 0.01950971
## GO:0004673 0.0005917603 0.03089794 0.02794020
## GO:0000155 0.0006221061 0.03089794 0.02794020
## geneID
## GO:0005887 sll0016/sll0094/sll0477/sll0478/sll0648/sll1124/sll1191/sll1192/sll1404/sll1405/sll1600/slr0210/slr0533/slr1147/slr1249/slr1317/slr1400/slr1647/slr1805/slr2045
## GO:0004673 sll0094/sll0790/sll1124/sll1888/sll1905/slr0210/slr0533/slr1147/slr1212/slr1285/slr1400/slr1731/slr1805
## GO:0000155 sll0094/sll0790/sll1124/sll1888/sll1905/slr0210/slr0487/slr0533/slr0829/slr1147/slr1212/slr1285/slr1400/slr1731/slr1805
## Count
## GO:0005887 20
## GO:0004673 13
## GO:0000155 15
## ID Description GeneRatio
## GO:0042651 GO:0042651 thylakoid membrane 36/371
## GO:0008137 GO:0008137 NADH dehydrogenase (ubiquinone) activity 13/371
## GO:0019684 GO:0019684 photosynthesis, light reaction 13/371
## GO:0048038 GO:0048038 quinone binding 13/371
## GO:0042773 GO:0042773 ATP synthesis coupled electron transport 10/371
## GO:0015986 GO:0015986 ATP synthesis coupled proton transport 8/371
## BgRatio pvalue p.adjust qvalue
## GO:0042651 94/2476 1.332696e-08 1.959064e-06 1.753548e-06
## GO:0008137 19/2476 1.822145e-07 1.339277e-05 1.198780e-05
## GO:0019684 20/2476 4.505429e-07 2.207660e-05 1.976065e-05
## GO:0048038 21/2476 1.023581e-06 3.761661e-05 3.367044e-05
## GO:0042773 14/2476 2.908849e-06 8.552015e-05 7.654865e-05
## GO:0015986 10/2476 8.100531e-06 1.984630e-04 1.776432e-04
## geneID
## GO:0042651 sll0519/sll0520/sll0522/sll0629/sll0689/sll1316/sll1322/sll1323/sll1324/sll1325/sll1326/sll1867/slr0228/slr0261/slr0782/slr0906/slr0927/slr1279/slr1280/slr1281/slr1291/slr1311/slr1329/slr1513/slr1623/slr1655/slr1796/sml0001/smr0006/smr0007/smr0008/smr0009/ssl2615/ssl3335/ssr1386/ssr3451
## GO:0008137 sll0027/sll0520/sll0521/sll0522/sll1732/sll1733/slr0261/slr1279/slr1280/slr1281/slr1291/slr2007/slr2009
## GO:0019684 sll0519/sll0520/sll0522/sll1867/slr0261/slr0906/slr0927/slr1279/slr1280/slr1281/slr1311/smr0006/ssr3451
## GO:0048038 sll0519/sll0520/sll0521/sll0522/sll1732/slr0261/slr1279/slr1280/slr1281/slr1623/slr2007/slr2009/ssr1386
## GO:0042773 sll0027/sll0522/sll1732/sll1733/slr1136/slr1291/slr2006/slr2007/slr2008/slr2009
## GO:0015986 sll1321/sll1322/sll1323/sll1324/sll1325/sll1326/slr1329/ssl2615
## Count
## GO:0042651 36
## GO:0008137 13
## GO:0019684 13
## GO:0048038 13
## GO:0042773 10
## GO:0015986 8go_functional_enrichment(CDS_result_dWT_dIF_0h, write=TRUE, path_up="output/enrichment/go_enrichment/dWT_dIF_0h_up.csv", path_down="output/enrichment/go_enrichment/dWT_dIF_0h_down.csv")## ID Description GeneRatio BgRatio pvalue p.adjust
## GO:0003677 GO:0003677 DNA binding 8/36 145/2423 0.0009785381 0.04892691
## qvalue
## GO:0003677 0.04223165
## geneID
## GO:0003677 sll0649/sll0856/sll1689/sll7106/slr0449/slr7049/slr7071/slr7092
## Count
## GO:0003677 8
## [1] ID Description GeneRatio BgRatio pvalue p.adjust
## [7] qvalue geneID Count
## <0 rows> (or 0-length row.names)go_functional_enrichment(CDS_result_dWT_dIF_1h, write=TRUE, path_up="output/enrichment/go_enrichment/dWT_dIF_1h_up.csv", path_down="output/enrichment/go_enrichment/dWT_dIF_1h_down.csv")## ID Description
## GO:0015979 GO:0015979 photosynthesis
## GO:0030089 GO:0030089 phycobilisome
## GO:0030096 GO:0030096 plasma membrane-derived thylakoid photosystem II
## GO:0042651 GO:0042651 thylakoid membrane
## GO:0009523 GO:0009523 photosystem II
## GO:0051607 GO:0051607 <NA>
## GeneRatio BgRatio pvalue p.adjust qvalue
## GO:0015979 28/211 85/2473 6.415827e-11 8.148100e-09 7.766527e-09
## GO:0030089 9/211 24/2473 8.504556e-05 5.400393e-03 5.147495e-03
## GO:0030096 9/211 27/2473 2.427170e-04 7.707963e-03 7.347002e-03
## GO:0042651 19/211 94/2473 2.427705e-04 7.707963e-03 7.347002e-03
## GO:0009523 6/211 13/2473 3.703654e-04 9.407282e-03 8.966742e-03
## GO:0051607 5/211 10/2473 7.616605e-04 1.612181e-02 1.536683e-02
## geneID
## GO:0015979 sll0047/sll0427/sll0629/sll1214/sll1398/sll1579/sll1580/slr0149/slr0506/slr1655/slr1739/slr1834/slr1835/slr1986/slr2051/slr2067/sml0001/sml0008/smr0004/smr0005/smr0006/smr0007/smr0008/ssl3093/ssr0330/ssr3304/ssr3383/ssr3451
## GO:0030089 sll1579/sll1580/slr0149/slr1986/slr2051/slr2067/slr7023/ssl3093/ssr3383
## GO:0030096 sll0427/sll1398/sll1867/slr0927/sml0001/smr0006/smr0007/smr0008/ssr3451
## GO:0042651 sll0047/sll0616/sll0629/sll1867/slr0927/slr1034/slr1655/slr1834/slr1835/slr1949/sml0001/sml0008/smr0004/smr0005/smr0006/smr0007/smr0008/ssl3335/ssr3451
## GO:0009523 sll0047/sll1867/slr0927/smr0006/smr0008/ssr3451
## GO:0051607 slr7016/slr7071/slr7092/ssr7017/ssr7093
## Count
## GO:0015979 28
## GO:0030089 9
## GO:0030096 9
## GO:0042651 19
## GO:0009523 6
## GO:0051607 5
## [1] ID Description GeneRatio BgRatio pvalue p.adjust
## [7] qvalue geneID Count
## <0 rows> (or 0-length row.names)kegg_functional_enrichment(CDS_result_dWT, write=TRUE, path_up="output/enrichment/kegg_enrichment/dWT_0h_1h_up.csv", path_down="output/enrichment/kegg_enrichment/dWT_0h_1h_down.csv")## Reading KEGG annotation online:
##
## Reading KEGG annotation online:## [1] ID Description GeneRatio BgRatio pvalue p.adjust
## [7] qvalue geneID Count
## <0 rows> (or 0-length row.names)
## ID Description GeneRatio BgRatio pvalue
## syn00190 syn00190 Oxidative phosphorylation 19/104 49/949 1.418619e-07
## p.adjust qvalue
## syn00190 6.525648e-06 6.421118e-06
## geneID
## syn00190 sll0223/sll0519/sll0520/sll0521/sll0522/sll1324/sll1325/sll1326/sll1327/sll1732/sll1733/slr0261/slr1279/slr1280/slr1281/slr1291/slr1330/slr1623/ssl2615
## Count
## syn00190 19kegg_functional_enrichment(CDS_result_dIF, write=TRUE, path_up="output/enrichment/kegg_enrichment/dIF_0h_1h_up.csv", path_down="output/enrichment/kegg_enrichment/dIF_0h_1h_down.csv")## [1] ID Description GeneRatio BgRatio pvalue p.adjust
## [7] qvalue geneID Count
## <0 rows> (or 0-length row.names)
## ID Description GeneRatio BgRatio pvalue
## syn00190 syn00190 Oxidative phosphorylation 23/172 49/949 1.760797e-06
## syn00195 syn00195 Photosynthesis 24/172 62/949 6.183843e-05
## p.adjust qvalue
## syn00190 9.860462e-05 9.638045e-05
## syn00195 1.731476e-03 1.692420e-03
## geneID
## syn00190 sll0027/sll0519/sll0520/sll0521/sll0522/sll1322/sll1323/sll1324/sll1325/sll1326/sll1732/sll1733/slr0261/slr1136/slr1137/slr1279/slr1280/slr1281/slr1291/slr1329/slr1623/ssl2615/ssr1386
## syn00195 sll0427/sll0629/sll1316/sll1322/sll1323/sll1324/sll1325/sll1326/sll1398/sll1796/sll1867/slr0906/slr0927/slr1311/slr1329/slr1655/slr1828/sml0001/smr0006/smr0007/smr0008/ssl0020/ssl2615/ssr3451
## Count
## syn00190 23
## syn00195 24kegg_functional_enrichment(CDS_result_dWT_dIF_0h, write=TRUE, path_up="output/enrichment/kegg_enrichment/dWT_dIF_0h_up.csv", path_down="output/enrichment/kegg_enrichment/dWT_dIF_0h_down.csv")## [1] ID Description GeneRatio BgRatio pvalue p.adjust
## [7] qvalue geneID Count
## <0 rows> (or 0-length row.names)## --> No gene can be mapped....## --> Expected input gene ID: sll1479,slr0301,slr1051,sll0745,slr0435,sll1213## --> return NULL...## NULLkegg_functional_enrichment(CDS_result_dWT_dIF_1h, write=TRUE, path_up="output/enrichment/kegg_enrichment/dWT_dIF_1h_up.csv", path_down="output/enrichment/kegg_enrichment/dWT_dIF_1h_down.csv")## ID Description GeneRatio BgRatio
## syn00195 syn00195 Photosynthesis 18/73 62/949
## syn00196 syn00196 Photosynthesis - antenna proteins 7/73 15/949
## pvalue p.adjust qvalue
## syn00195 1.485540e-07 6.684929e-06 6.684929e-06
## syn00196 4.699591e-05 1.057408e-03 1.057408e-03
## geneID
## syn00195 sll0427/sll0629/sll1182/sll1398/sll1867/slr0927/slr1655/slr1739/slr1834/slr1835/sml0001/sml0008/smr0004/smr0005/smr0006/smr0007/smr0008/ssr3451
## syn00196 sll1579/sll1580/slr1986/slr2051/slr2067/ssl3093/ssr3383
## Count
## syn00195 18
## syn00196 7
## [1] ID Description GeneRatio BgRatio pvalue p.adjust
## [7] qvalue geneID Count
## <0 rows> (or 0-length row.names)3.3 Gene Set Enrichment Analyses
dWT_0h_1h_go_gsea <- go_gsea(CDS_result_dWT, write=TRUE, path="output/enrichment/go_gsea/dWT_0h_1h_go_gsea.csv")## preparing geneSet collections...## GSEA analysis...## Warning in fgseaMultilevel(...): For some pathways, in reality P-values are less
## than 1e-10. You can set the `eps` argument to zero for better estimation.## leading edge analysis...## done...## ID Description setSize
## GO:0008137 GO:0008137 NADH dehydrogenase (ubiquinone) activity 19
## GO:0042773 GO:0042773 ATP synthesis coupled electron transport 14
## GO:0048038 GO:0048038 quinone binding 21
## GO:0055114 GO:0055114 oxidation-reduction process 201
## GO:0019684 GO:0019684 photosynthesis, light reaction 20
## GO:0015977 GO:0015977 carbon fixation 13
## enrichmentScore NES pvalue p.adjust qvalues
## GO:0008137 -0.9218843 -2.682260 1.000000e-10 1.760000e-08 1.410526e-08
## GO:0042773 -0.9163969 -2.464829 2.748115e-09 2.418341e-07 1.938144e-07
## GO:0048038 -0.8397445 -2.465698 8.177258e-09 4.797325e-07 3.844746e-07
## GO:0055114 -0.4423196 -1.961674 3.030759e-07 1.333534e-05 1.068741e-05
## GO:0019684 -0.7784678 -2.261980 4.913748e-06 1.729639e-04 1.386194e-04
## GO:0015977 -0.8257245 -2.165586 1.894119e-05 5.556084e-04 4.452842e-04
## rank leading_edge
## GO:0008137 120 tags=68%, list=2%, signal=67%
## GO:0042773 98 tags=64%, list=2%, signal=63%
## GO:0048038 120 tags=62%, list=2%, signal=61%
## GO:0055114 899 tags=28%, list=18%, signal=24%
## GO:0019684 120 tags=40%, list=2%, signal=39%
## GO:0015977 266 tags=62%, list=5%, signal=58%dIF_0h_1h_go_gsea <- go_gsea(CDS_result_dIF, write=TRUE, path="output/enrichment/go_gsea/dIF_0h_1h_go_gsea.csv")## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...## ID Description setSize
## GO:0055114 GO:0055114 oxidation-reduction process 202
## GO:0042651 GO:0042651 thylakoid membrane 94
## GO:0015979 GO:0015979 photosynthesis 85
## GO:0006412 GO:0006412 translation 89
## GO:0003735 GO:0003735 structural constituent of ribosome 55
## GO:0005840 GO:0005840 ribosome 54
## enrichmentScore NES pvalue p.adjust qvalues
## GO:0055114 -0.4276135 -2.112589 6.073483e-10 6.040886e-08 4.443953e-08
## GO:0042651 -0.5499966 -2.399096 6.864643e-10 6.040886e-08 4.443953e-08
## GO:0015979 -0.5494639 -2.369990 8.059436e-09 4.084434e-07 3.004697e-07
## GO:0006412 -0.5350417 -2.321256 9.282804e-09 4.084434e-07 3.004697e-07
## GO:0003735 -0.6065864 -2.385362 5.033700e-08 1.771863e-06 1.303463e-06
## GO:0005840 -0.6003029 -2.358980 1.347555e-07 3.952828e-06 2.907882e-06
## rank leading_edge
## GO:0055114 1403 tags=42%, list=25%, signal=33%
## GO:0042651 1276 tags=59%, list=23%, signal=46%
## GO:0015979 1601 tags=61%, list=29%, signal=44%
## GO:0006412 1724 tags=71%, list=31%, signal=49%
## GO:0003735 1667 tags=80%, list=30%, signal=56%
## GO:0005840 1667 tags=80%, list=30%, signal=56%dWT_dIF_0h_go_gsea <- go_gsea(CDS_result_dWT_dIF_0h, write=TRUE, path="output/enrichment/go_gsea/dWT_dIF_0h_go_gsea.csv")## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...## ID Description setSize
## GO:0004519 GO:0004519 endonuclease activity 24
## GO:0006355 GO:0006355 regulation of transcription, DNA-templated 135
## GO:0003677 GO:0003677 DNA binding 145
## enrichmentScore NES pvalue p.adjust qvalues rank
## GO:0004519 0.7330815 2.018810 3.323554e-05 0.005782985 0.005527596 450
## GO:0006355 0.4386987 1.691613 4.164360e-04 0.028494167 0.027235804 1053
## GO:0003677 0.4361410 1.687431 4.912787e-04 0.028494167 0.027235804 739
## leading_edge
## GO:0004519 tags=46%, list=10%, signal=41%
## GO:0006355 tags=36%, list=24%, signal=28%
## GO:0003677 tags=27%, list=17%, signal=23%dWT_dIF_1h_go_gsea <- go_gsea(CDS_result_dWT_dIF_1h, write=TRUE, path="output/enrichment/go_gsea/dWT_dIF_1h_go_gsea.csv")## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...## ID Description setSize
## GO:0015979 GO:0015979 photosynthesis 85
## GO:0042651 GO:0042651 thylakoid membrane 94
## GO:0006412 GO:0006412 translation 89
## GO:0003735 GO:0003735 structural constituent of ribosome 55
## GO:0005840 GO:0005840 ribosome 54
## GO:0030089 GO:0030089 phycobilisome 24
## enrichmentScore NES pvalue p.adjust qvalues rank
## GO:0015979 0.5607895 2.468703 3.516686e-10 6.189368e-08 5.441610e-08 864
## GO:0042651 0.4743434 2.098625 4.106142e-07 3.613405e-05 3.176857e-05 1375
## GO:0006412 0.4346562 1.920062 1.798293e-05 1.054999e-03 9.275406e-04 1634
## GO:0003735 0.4972601 1.983787 4.788636e-05 2.107000e-03 1.852446e-03 1563
## GO:0005840 0.4963100 1.980938 7.859357e-05 2.590323e-03 2.277378e-03 1563
## GO:0030089 0.6440629 2.125906 8.830648e-05 2.590323e-03 2.277378e-03 995
## leading_edge
## GO:0015979 tags=47%, list=16%, signal=40%
## GO:0042651 tags=46%, list=26%, signal=34%
## GO:0006412 tags=57%, list=31%, signal=40%
## GO:0003735 tags=65%, list=30%, signal=47%
## GO:0005840 tags=65%, list=30%, signal=46%
## GO:0030089 tags=58%, list=19%, signal=48%dWT_0h_1h_kegg_gsea <- kegg_gsea(CDS_result_dWT, write=TRUE, path="output/enrichment/kegg_gsea/dWT_0h_1h_kegg_gsea.csv")## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...## ID Description setSize enrichmentScore NES
## syn00190 syn00190 Oxidative phosphorylation 49 -0.7011069 -2.46859
## pvalue p.adjust qvalues rank
## syn00190 8.115011e-09 4.950156e-07 4.014795e-07 896
## leading_edge
## syn00190 tags=59%, list=18%, signal=49%dIF_0h_1h_kegg_gsea <- kegg_gsea(CDS_result_dIF, write=TRUE, path="output/enrichment/kegg_gsea/dIF_0h_1h_kegg_gsea.csv")## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...## ID Description setSize enrichmentScore
## syn00190 syn00190 Oxidative phosphorylation 49 -0.6817810
## syn00195 syn00195 Photosynthesis 62 -0.6048131
## syn03010 syn03010 Ribosome 53 -0.6033970
## syn00240 syn00240 Pyrimidine metabolism 25 -0.6715112
## syn01110 syn01110 Biosynthesis of secondary metabolites 285 -0.2879301
## syn03070 syn03070 Bacterial secretion system 14 -0.6774017
## NES pvalue p.adjust qvalues rank
## syn00190 -2.664246 2.962910e-10 1.807375e-08 1.341106e-08 1115
## syn00195 -2.495725 6.553581e-09 1.998842e-07 1.483179e-07 1608
## syn03010 -2.381502 8.731977e-08 1.775502e-06 1.317456e-06 1667
## syn00240 -2.213964 3.859505e-05 5.885745e-04 4.367335e-04 985
## syn01110 -1.483742 3.723699e-04 4.542913e-03 3.370928e-03 1417
## syn03070 -1.907194 2.445821e-03 2.486585e-02 1.845093e-02 1788
## leading_edge
## syn00190 tags=63%, list=20%, signal=51%
## syn00195 tags=73%, list=29%, signal=52%
## syn03010 tags=79%, list=30%, signal=56%
## syn00240 tags=52%, list=18%, signal=43%
## syn01110 tags=33%, list=26%, signal=26%
## syn03070 tags=100%, list=32%, signal=68%dWT_dIF_0h_kegg_gsea <- kegg_gsea(CDS_result_dWT_dIF_0h, write=TRUE, path="output/enrichment/kegg_gsea/dWT_dIF_0h_kegg_gsea.csv")## preparing geneSet collections...## GSEA analysis...## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"dWT_dIF_1h_kegg_gsea <- kegg_gsea(CDS_result_dWT_dIF_1h, write=TRUE, path="output/enrichment/kegg_gsea/dWT_dIF_1h_kegg_gsea.csv")## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...## ID Description setSize enrichmentScore
## syn00195 syn00195 Photosynthesis 62 0.5836453
## syn03010 syn03010 Ribosome 53 0.5039717
## syn00196 syn00196 Photosynthesis - antenna proteins 15 0.6968561
## syn00790 syn00790 Folate biosynthesis 17 -0.6603830
## NES pvalue p.adjust qvalues rank
## syn00195 2.367202 2.010784e-08 1.226578e-06 1.142972e-06 1044
## syn03010 1.971284 1.157646e-04 3.530821e-03 3.290152e-03 1563
## syn00196 2.037089 2.341332e-04 4.760709e-03 4.436208e-03 948
## syn00790 -1.835465 1.937164e-03 2.954175e-02 2.752812e-02 970
## leading_edge
## syn00195 tags=53%, list=20%, signal=43%
## syn03010 tags=68%, list=30%, signal=48%
## syn00196 tags=73%, list=18%, signal=60%
## syn00790 tags=47%, list=18%, signal=39%3.4 Plots of GSEA
gseaplot2(dWT_0h_1h_go_gsea, geneSetID=1:6)
p <- gseaplot2(dWT_0h_1h_kegg_gsea, geneSetID =1)
p
ggsave("output/enrichment/Plots/dWT_0h_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")p <- gseaplot2(dIF_0h_1h_kegg_gsea, geneSetID =1:6)
p
ggsave("output/enrichment/Plots/dIF_0h_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn00190", 1) # https://www.kegg.jp/kegg-bin/show_pathway?syn00190/sll0026/slr2083/sll0223/sll1899/sll1262/sll1327/slr1622/sll1223/slr1279/slr1329/slr1137/sll0027/sll1326/ssr1386/sll0519/slr0261/sll1322/slr1136/slr1623/slr1280/sll1324/ssl2615/sll1325/sll1323/sll0522/sll0520/slr1281/sll0521/sll1733/slr1291/sll1732
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn00195", 2) # https://www.kegg.jp/kegg-bin/show_pathway?syn00195/sll1382/ssr0390/slr0343/sll0849/slr0342/sml0005/ssl2598/smr0005/slr0737/slr1645/slr1835/slr1643/smr0003/sll0819/sll1182/sll1327/slr1185/slr1834/smr0010/smr0004/sll1317/sll0427/smr0008/sll1796/slr1329/sll1316/slr1655/slr0906/sml0001/sll1326/slr1311/smr0006/smr0007/ssr3451/sll0629/slr1828/slr0927/sll1322/sll1324/ssl2615/sll1325/sll1323/sll1867/ssl0020/sll1398
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn03010", 3) # https://www.kegg.jp/kegg-bin/show_pathway?syn03010/slr1356/slr1678/sll1800/sll1260/sll1822/sll1801/sll1096/sll1821/sll1799/sll1808/ssr1604/sll1824/sll1816/sll1802/ssr1398/sll1804/sll1813/sll1806/ssl1784/sll1097/sll1807/ssr1399/ssl3437/sll1809/ssl3432/ssr0482/sll1746/sll1812/sll1817/ssl3436/ssr1736/sll1803/sll1805/sll1811/sll1244/sll0767/ssl1426/sll1743/sll1740/sll1745/slr1984/smr0011
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn00240", 4) # https://www.kegg.jp/kegg-bin/show_pathway?syn00240/sll1459/sll1108/slr1476/sll1635/slr1616/sll0144/slr1164/sll1258/sll0368/sll1035/sll1852/sll0370/slr0591
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn01110", 5) # https://www.kegg.jp/kegg-bin/show_pathway?syn01110/sll0902/sll0329/sll1558/slr1022/slr0738/slr0186/slr1289/slr1898/slr0049/slr0506/slr1887/slr1019/slr1808/slr0886/sll1444/slr1790/slr1867/slr0879/slr1055/sll1127/sll1468/sll1185/slr1226/sll1538/sll1172/sll1688/slr2136/slr0444/slr0293/slr1737/sll0513/sll1994/sll1899/sll0179/sll0901/sll1459/sll1108/sll0033/slr2072/sll1479/slr1349/slr0585/slr2130/sll1056/sll1747/slr0526/slr1124/sll1349/sll0158/sll1908/slr1176/sll1363/sll1091/slr0018/slr1510/slr1096/sll0807/sll0421/sll0480/slr1165/sll1342/slr0657/slr1743/slr0611/slr0985/slr1934/slr0966/slr0940/slr0116/slr1736/sll1815/slr1842/sll1797/slr1517/slr2094/slr0394/slr1945/sll0373/slr1843/sll0735/slr2023/sll0418/slr1030/sll0422/sll0018/sll1214/slr0749/slr0348/slr0012/slr0009/sll1852/sll0927/slr0493/sll0019/ssl2084
browseKEGGNew_3(dIF_0h_1h_kegg_gsea, "syn03070", 6) # https://www.kegg.jp/kegg-bin/show_pathway?syn03070/slr1471/sll1180/ssr3307/sll1814/sll1181/slr1531/sll0194/slr1046/sll0616/slr0774/slr0775/ssl2823/ssl3335
browseKEGGNew_3(dWT_0h_1h_kegg_gsea, "syn00190", 1) # https://www.kegg.jp/kegg-bin/show_pathway?syn00190/sll1899/slr2082/sll0026/sll1262/sll1322/slr1622/sll1484/slr1136/sll1323/slr1137/ssl2615/slr1330/slr1623/sll1327/sll1326/sll1324/sll1325/slr0261/slr1279/sll0223/sll0519/slr1280/sll0521/sll0520/sll0522/slr1281/slr1291/sll1733/sll1732p <- gseaplot2(dWT_dIF_1h_kegg_gsea, geneSetID=1:4)
p
ggsave("output/enrichment/Plots/dWT_dIF_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")
p <- gseaplot2(dWT_dIF_1h_go_gsea, geneSetID=c(2,8,3))
p
ggsave("output/enrichment/Plots/dWT_dIF_1h_GO.pdf", plot=p, width=15, height=12, units="cm")
p <- gseaplot2(dWT_dIF_1h_kegg_gsea, geneSetID=c(1,3,2,4))
p
ggsave("output/enrichment/Plots/dWT_dIF_1h_KEGG.pdf", plot=p, width=15, height=12, units="cm")browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn00195", 1) # https://www.kegg.jp/kegg-bin/show_pathway?syn00195/smr0005/slr1655/smr0004/sll0629/slr1834/slr1835/ssr3451/sll1182/smr0007/smr0006/smr0008/slr1739/sml0008/sll1867/sll0427/slr0927/sml0001/sll1398/sll1323/sll1382/slr1311/ssl0020/smr0010/sml0002/slr0737/sll0819/sll1322/slr0150/ssl0563/sll1325/slr0906/sll0849/sll1324
browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn03010", 2) # https://www.kegg.jp/kegg-bin/show_pathway?syn03010/sll1740/smr0011/slr1984/ssr0482/sll1811/sll1244/sll1743/ssl1784/ssr1604/sll1803/sll1812/ssl1426/sll1824/sll1813/sll1806/sll1807/ssr1736/ssl3437/sll1096/sll1805/ssl3432/sll1809/sll1810/ssl3436/sll1802/ssr1398/sll1804/slr1678/sll1799/ssr1399/sll1800/sll1808/sll1821/sll1097/sll1801/sll1260
browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn00196", 3) # https://www.kegg.jp/kegg-bin/show_pathway?syn00196/slr1986/slr2067/sll1579/sll1580/ssr3383/ssl3093/slr2051/sll1577/sll1578/slr0335/sll0928
browseKEGGNew_3(dWT_dIF_1h_kegg_gsea, "syn00790", 4) # https://www.kegg.jp/kegg-bin/show_pathway?syn00790/slr0527/slr1093/slr0078/slr0887/slr0900/slr0902/slr0901/slr09033.5 GSEA for RNA features
First, a data frame in which feature names are assigned to their feature type (CDS, ncRNA, …) has to be created and the respective info read in.
df_featureType <- as.data.frame(cbind(feature_type=as.character(features$type), feature_name=features$locus_tag))
feature_type_gsea <- function(DESEq2_dataframe, write=FALSE, path=""){
locus_tags <- row.names(DESEq2_dataframe)
geneList <- DESEq2_dataframe$log2FoldChange
names(geneList) <- locus_tags
geneList = sort(geneList, decreasing = TRUE)
set.seed(42)
go_gsea_object <- GSEA(geneList, TERM2GENE = df_featureType, seed=TRUE)
print(head(go_gsea_object)[,1:10])
if(write){
write.csv(go_gsea_object, path)
}
return(go_gsea_object)
}
features_gsea_0h <- feature_type_gsea(CDS_result_dWT_dIF_0h, TRUE, "output/enrichment/GSEA_RNAfeatures_dWT-dIF_0h.csv")## preparing geneSet collections...## GSEA analysis...## Warning in fgseaMultilevel(...): For some pathways, in reality P-values are less
## than 1e-10. You can set the `eps` argument to zero for better estimation.## leading edge analysis...## done...## ID Description setSize enrichmentScore NES pvalue
## 5UTR 5UTR 5UTR 272 -0.4118204 -2.119681 1.000000e-10
## tRNA tRNA tRNA 35 0.7360567 2.283511 5.012588e-07
## ncRNA ncRNA ncRNA 223 -0.2909000 -1.443708 2.227298e-03
## 3UTR 3UTR 3UTR 25 0.6002223 1.745930 6.680201e-03
## p.adjust qvalues rank leading_edge
## 5UTR 4.000000e-10 NA 846 tags=41%, list=20%, signal=35%
## tRNA 1.002518e-06 NA 652 tags=63%, list=15%, signal=54%
## ncRNA 2.969731e-03 NA 627 tags=28%, list=15%, signal=25%
## 3UTR 6.680201e-03 NA 363 tags=40%, list=8%, signal=37%features_gsea_1h <- feature_type_gsea(CDS_result_dWT_dIF_1h, TRUE, "output/enrichment/GSEA_RNAfeatures_dWT-dIF_1h.csv")## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...## ID Description setSize enrichmentScore NES pvalue
## tRNA tRNA tRNA 38 0.6716029 2.459963 2.409659e-08
## 3UTR 3UTR 3UTR 29 0.5561780 1.910522 7.894904e-04
## p.adjust qvalues rank leading_edge
## tRNA 7.228978e-08 2.536484e-08 489 tags=53%, list=9%, signal=48%
## 3UTR 1.184236e-03 4.155212e-04 1062 tags=48%, list=20%, signal=39%p <- gseaplot2(features_gsea_0h, geneSetID=1:4)
p
p <- gseaplot2(features_gsea_1h, geneSetID=1:2)
p
3.6 GSEA for base means
go_gsea_baseMeans <- function(DESEq2_dataframe, write=FALSE, path=""){
locus_tags <- row.names(DESEq2_dataframe)
geneList <- DESEq2_dataframe$baseMean
names(geneList) <- locus_tags
geneList = sort(geneList, decreasing = TRUE)
set.seed(42)
go_gsea_object <- GSEA(geneList, TERM2GENE = term_to_gene, TERM2NAME=term_to_name, seed=TRUE)
tryCatch({
print(head(go_gsea_object)[,1:10])
if(write){
write.csv(go_gsea_object, path)
}
return(go_gsea_object)}, error=function(e){
print("nothing enriched")
})
}
kegg_gsea_baseMean <- function(DESEq2_dataframe, write=FALSE, path=""){
locus_tags <- row.names(DESEq2_dataframe)
geneList <- DESEq2_dataframe$baseMean
names(geneList) <- locus_tags
geneList = sort(geneList, decreasing = TRUE)
set.seed(42)
kegg_gsea_object <- gseKEGG(geneList, organism="syn", minGSSize=10, pvalueCutoff = 0.05, seed=TRUE)
tryCatch({
print(head(kegg_gsea_object)[,1:10])
if(write){
write.csv(kegg_gsea_object, path)
}
return(kegg_gsea_object)}, error=function(e){
print("nothing enriched")
})
}3.7 Without taking width into account
When width of features is not taken into account, there is no specific enrichment of certain terms in more highly expressed features. One exception is the GO term “transposase activity”, which is depleted when analyzing the comparison of dWT and dIF at time point 0h.
go_gsea_baseMeans(CDS_result_dIF)## preparing geneSet collections...## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".## Warning in fgseaMultilevel(...): There were 2 pathways for which P-values were
## not calculated properly due to unbalanced (positive and negative) gene-level
## statistic values. For such pathways pval, padj, NES, log2err are set to NA. You
## can try to increase the value of the argument nPermSimple (for example set it
## nPermSimple = 10000)## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"go_gsea_baseMeans(CDS_result_dWT)## preparing geneSet collections...
## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There were 2 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple =
## 10000)## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"go_gsea_baseMeans(CDS_result_dWT_dIF_0h)## preparing geneSet collections...
## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".## leading edge analysis...## done...## ID Description setSize enrichmentScore NES
## GO:0004803 GO:0004803 transposase activity 11 -0.7504063 -2.409325
## pvalue p.adjust qvalues rank
## GO:0004803 1.316033e-05 0.002289897 0.002271888 1087
## leading_edge
## GO:0004803 tags=100%, list=25%, signal=75%## #
## # Gene Set Enrichment Analysis
## #
## #...@organism UNKNOWN
## #...@setType UNKNOWN
## #...@geneList Named num [1:4318] 209545 165602 81721 65565 59371 ...
## - attr(*, "names")= chr [1:4318] "ncl1680" "ncr0700" "ncr0075" "sll1578" ...
## #...nPerm
## #...pvalues adjusted by 'BH' with cutoff <0.05
## #...1 enriched terms found
## 'data.frame': 1 obs. of 11 variables:
## $ ID : chr "GO:0004803"
## $ Description : chr "transposase activity"
## $ setSize : int 11
## $ enrichmentScore: num -0.75
## $ NES : num -2.41
## $ pvalue : num 1.32e-05
## $ p.adjust : num 0.00229
## $ qvalues : num 0.00227
## $ rank : num 1087
## $ leading_edge : chr "tags=100%, list=25%, signal=75%"
## $ core_enrichment: chr "sll1861/sll1710/ssl3649/sll1716/slr0800/slr2096/slr1903/slr5040/sll1792/sll1157"
## #...Citation
## Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He.
## clusterProfiler: an R package for comparing biological themes among
## gene clusters. OMICS: A Journal of Integrative Biology
## 2012, 16(5):284-287go_gsea_baseMeans(CDS_result_dWT_dIF_1h)## preparing geneSet collections...## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".
## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : There were 2 pathways for which P-values were not calculated
## properly due to unbalanced (positive and negative) gene-level statistic values.
## For such pathways pval, padj, NES, log2err are set to NA. You can try to
## increase the value of the argument nPermSimple (for example set it nPermSimple =
## 10000)## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"kegg_gsea_baseMean(CDS_result_dIF)## preparing geneSet collections...
## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"kegg_gsea_baseMean(CDS_result_dWT)## preparing geneSet collections...
## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"kegg_gsea_baseMean(CDS_result_dWT_dIF_0h)## preparing geneSet collections...
## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"kegg_gsea_baseMean(CDS_result_dWT_dIF_1h)## preparing geneSet collections...
## GSEA analysis...## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize,
## gseaParam, : All values in the stats vector are greater than zero and scoreType
## is "std", maybe you should switch to scoreType = "pos".## no term enriched under specific pvalueCutoff...## [1] "nothing enriched"3.8 When taking width of features into account
norm_counts <- counts(ddsMat_CDS, normalized=TRUE)
mean_norm_counts <- apply(norm_counts, 1, mean)
mean_norm_counts_CDS <- subset(mean_norm_counts, names(mean_norm_counts) %in% CDS_features$locus_tag)
CDS_feat_df <- as.data.frame(CDS_features)
row.names(CDS_feat_df) <- CDS_feat_df$locus_tag
mean_norm_counts_CDS_width <- mean_norm_counts_CDS/CDS_feat_df[names(mean_norm_counts_CDS),]$widthmean_norm_counts_CDS_width <- sort(mean_norm_counts_CDS_width, decreasing=TRUE)
set.seed(42)
go_gsea_baseMeans_width <- GSEA(mean_norm_counts_CDS_width, TERM2GENE=term_to_gene, TERM2NAME = term_to_name, nPermSimple=10000, seed=TRUE)## preparing geneSet collections...## GSEA analysis...## Warning in fgseaMultilevel(...): There were 3 pathways for which P-values were
## not calculated properly due to unbalanced (positive and negative) gene-level
## statistic values. For such pathways pval, padj, NES, log2err are set to NA. You
## can try to increase the value of the argument nPermSimple (for example set it
## nPermSimple = 100000)## leading edge analysis...## done...set.seed(42)
kegg_gsea_baseMeans_width <- gseKEGG(mean_norm_counts_CDS_width, organism="syn", minGSSize=10, pvalueCutoff = 0.05, nPermSimple=10000, seed=TRUE)## preparing geneSet collections...## GSEA analysis...## leading edge analysis...## done...go_gsea_baseMeans_width[,1:10]## ID Description setSize enrichmentScore
## GO:0015979 GO:0015979 photosynthesis 86 0.9254062
## GO:0042651 GO:0042651 thylakoid membrane 96 0.8945113
## GO:0018298 GO:0018298 protein-chromophore linkage 31 0.9451642
## GO:0008272 GO:0008272 sulfate transport 12 -0.5812176
## NES pvalue p.adjust qvalues rank
## GO:0015979 1.565034 4.810606e-07 8.466667e-05 8.102074e-05 189
## GO:0042651 1.509690 1.417049e-04 1.247003e-02 1.193305e-02 196
## GO:0018298 1.621643 8.657711e-04 4.930474e-02 4.718157e-02 127
## GO:0008272 -1.811787 1.120562e-03 4.930474e-02 4.718157e-02 1547
## leading_edge
## GO:0015979 tags=56%, list=5%, signal=54%
## GO:0042651 tags=50%, list=5%, signal=49%
## GO:0018298 tags=45%, list=3%, signal=44%
## GO:0008272 tags=100%, list=42%, signal=58%write.csv(go_gsea_baseMeans_width, "output/enrichment/GO_GSEA_baseMeans_perWidth.csv")
kegg_gsea_baseMeans_width[,1:10]## ID Description setSize enrichmentScore NES pvalue
## syn00195 syn00195 Photosynthesis 63 0.9294204 1.580491 2.778544e-05
## p.adjust qvalues rank leading_edge
## syn00195 0.001722697 0.001722697 195 tags=70%, list=5%, signal=67%write.csv(kegg_gsea_baseMeans_width, "output/enrichment/KEGG_GSEA_baseMeans_perWidth.csv")p <- gseaplot2(go_gsea_baseMeans_width, geneSetID=c(1,2,3,4))
p
ggsave("output/enrichment/Plots/GO_GSEA_baseMeans_perWidth.pdf", plot=p, width=15, height=12, units="cm")p <- gseaplot2(kegg_gsea_baseMeans_width, geneSetID=c(1))
p
ggsave("output/enrichment/Plots/KEGG_GSEA_baseMeans_perWidth.pdf", plot=p, width=15, height=12, units="cm")Session info
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 11 (bullseye)
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=de_DE.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=de_DE.UTF-8 LC_COLLATE=de_DE.UTF-8
## [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=de_DE.UTF-8
## [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] forcats_0.5.1 stringr_1.4.0
## [3] dplyr_1.0.5 purrr_0.3.4
## [5] readr_1.4.0 tidyr_1.1.3
## [7] tibble_3.1.0 tidyverse_1.3.0
## [9] rtracklayer_1.50.0 enrichplot_1.10.2
## [11] clusterProfiler_3.18.1 RColorBrewer_1.1-2
## [13] pheatmap_1.0.12 ggpubr_0.4.0
## [15] ggrepel_0.9.1 ggplot2_3.3.3
## [17] DESeq2_1.30.1 SummarizedExperiment_1.20.0
## [19] Biobase_2.50.0 MatrixGenerics_1.2.1
## [21] matrixStats_0.58.0 GenomicRanges_1.42.0
## [23] GenomeInfoDb_1.26.2 IRanges_2.24.1
## [25] S4Vectors_0.28.1 BiocGenerics_0.36.0
## [27] knitr_1.31
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 shadowtext_0.0.7 backports_1.2.1
## [4] fastmatch_1.1-0 plyr_1.8.6 igraph_1.2.6
## [7] splines_4.0.4 BiocParallel_1.24.1 digest_0.6.27
## [10] htmltools_0.5.1.1 GOSemSim_2.16.1 viridis_0.5.1
## [13] GO.db_3.12.1 fansi_0.4.2 magrittr_2.0.1
## [16] memoise_2.0.0 openxlsx_4.2.3 Biostrings_2.58.0
## [19] annotate_1.68.0 graphlayouts_0.7.1 modelr_0.1.8
## [22] colorspace_2.0-0 rvest_1.0.0 blob_1.2.1
## [25] haven_2.3.1 xfun_0.22 crayon_1.4.1
## [28] RCurl_1.98-1.2 jsonlite_1.7.2 scatterpie_0.1.5
## [31] genefilter_1.72.1 survival_3.2-7 glue_1.4.2
## [34] polyclip_1.10-0 gtable_0.3.0 zlibbioc_1.36.0
## [37] XVector_0.30.0 DelayedArray_0.16.2 car_3.0-10
## [40] abind_1.4-5 scales_1.1.1 DOSE_3.16.0
## [43] DBI_1.1.1 rstatix_0.7.0 Rcpp_1.0.6
## [46] viridisLite_0.3.0 xtable_1.8-4 foreign_0.8-81
## [49] bit_4.0.4 httr_1.4.2 fgsea_1.16.0
## [52] ellipsis_0.3.1 pkgconfig_2.0.3 XML_3.99-0.5
## [55] farver_2.1.0 sass_0.3.1 dbplyr_2.1.0
## [58] locfit_1.5-9.4 utf8_1.1.4 labeling_0.4.2
## [61] tidyselect_1.1.0 rlang_0.4.10 reshape2_1.4.4
## [64] AnnotationDbi_1.52.0 munsell_0.5.0 cellranger_1.1.0
## [67] tools_4.0.4 cachem_1.0.4 cli_2.3.1
## [70] downloader_0.4 generics_0.1.0 RSQLite_2.2.3
## [73] broom_0.7.5 evaluate_0.14 fastmap_1.1.0
## [76] yaml_2.2.1 bit64_4.0.5 fs_1.5.0
## [79] tidygraph_1.2.0 zip_2.1.1 ggraph_2.0.5
## [82] xml2_1.3.2 DO.db_2.9 rstudioapi_0.13
## [85] compiler_4.0.4 curl_4.3 ggsignif_0.6.1
## [88] reprex_1.0.0 tweenr_1.0.1 geneplotter_1.68.0
## [91] bslib_0.2.4 stringi_1.5.3 highr_0.8
## [94] lattice_0.20-41 Matrix_1.3-2 vctrs_0.3.6
## [97] pillar_1.5.1 lifecycle_1.0.0 BiocManager_1.30.10
## [100] jquerylib_0.1.3 data.table_1.14.0 cowplot_1.1.1
## [103] bitops_1.0-6 qvalue_2.22.0 R6_2.5.0
## [106] gridExtra_2.3 rio_0.5.26 MASS_7.3-53.1
## [109] assertthat_0.2.1 withr_2.4.1 GenomicAlignments_1.26.0
## [112] Rsamtools_2.6.0 GenomeInfoDbData_1.2.4 hms_1.0.0
## [115] grid_4.0.4 rmarkdown_2.7 rvcheck_0.1.8
## [118] carData_3.0-4 ggforce_0.3.3 lubridate_1.7.10